#----this script define the various functions used for logistic model-----
# test on 6/10/2013

library(maxLik)

#-------This function return the P of the 3-parameter logistics model------

P_3pl<-function(theta,a,b,c){
  	DL=1.7*a*(theta-b);
        P=c+(1-c)*(exp(DL)/(1+exp(DL)));
	return(P);
}


#---------------------------------------------------------------------------
#info_3pl is the Information function for three parameter logistic model.
#	INFO(t,a,b,c) provides a vector of information values corresponding to the 
#	theta values in the t-vector.  The item parameters for the item are given by 
#	a, b, and c.  The model uses D=1.7 to maximize the similarity to the normal ogive
#	model.
#----------------------------------------------------------------------------------

info_3pl<-function(theta,a,b,c){

        DL=1.7*a*(theta-b);
        P=c+(1-c)*(exp(DL)/(1+exp(DL)));
        Q=1-P;
        i=1.7*1.7*a*a*(Q/P)*((P-c)/(1-c))^2;
}


#-----------------------------------------------------------------------------
#	testinfo_3pl define the Information function for three parameter logistic model.
#	INFO(t,parm) provides a vector of information values corresponding to the 
#	theta values in the t-vector.  The item parameters for the item are given by 
#	parm.  The model uses D=1.7 to maximize the similarity to the normal ogive
#	model.
#------------------------------------------------------------------------------------

testinfo_3pl<-function(theta,a,b,c){
          k=length(a);
	  m=length(theta);
          inf=matrix(0.,nrow=k,ncol=m)       
	  for(i in 1:k){
          aa=a[i];
          bb=b[i];
          cc=c[i];         
          for(j in 1:m){
            DL=1.7*aa*(theta[j]-bb);
            P=cc+(1-cc)*(exp(DL)/(1+exp(DL)));
            Q=1-P;
            inf[i,j]=1.7*1.7*aa*aa*(Q/P)*((P-cc)/(1-cc))^2;
         }
        }

        tinfo=colSums(inf)
}


#---------------------------------------------------------
#calculate the mean P of 3 parameter logistic model      -
#input is the theta, frequency of each theta, a, b and c - 
#---------------------------------------------------------

P_mean_MC<-function(theta,freq,a,b,c){
        s=0.;
        for(i in 1:length(theta)){
    		s=s+P_3pl(theta[i],a,b,c)*freq[i]
	} 
	res = s/sum(freq)
}


#-------------------------------------------------
# calculate the TCC for 3 parameter logistic model
#-------------------------------------------------

tcc_3pl<-function(theta,a,b,c){
   Ntheta=length(theta)
   tcc=matrix(0,nrow=Ntheta)
   for(i in 1:Ntheta){
	tcc[i]=sum(P_3pl(theta[i],a,b,c))
   }
   return(tcc);
}

#-----------------------------------------------------------------
#standard error of 3pl
#-----------------------------------------------------------------
sderr_3pl<-function(theta,a,b,c){
      res=1./sqrt(testinfo_3pl(theta,a,b,c));
}


#------------------------------------------------------------
#P_mean_MC from integration, by assumine theta as N(mu,sigma^2)
#-----------------------------------------------------------

pmean <- function(a,b,c,mu,sigma){
            pmean_int <- function(theta,a,b,c,mu,sigma){
                 res=P_3pl(a,b,c,theta)*exp(-(theta-mu)^2/2./sigma^2)/sqrt(2*pi)/sigma;
                 return(res);
               }

            res=integrate(pmean_int,lower=-Inf, upper=Inf,a=a,b=b,c=c,mu=mu,sigma=sigma,rel.tol = .Machine$double.eps^0.5)$val
            return(res);
          }

pmean_diff <- function(a,b,c){
      mu1=1;
      mu2=0;
      sigma1=0.3;
      sigma2=0.3;
      #a <- param[1]
      #b <- param[2]
      #c <- param[3]  
      res=pmean(a,b,c,mu1,sigma1)-pmean(a,b,c,mu2,sigma2);
      return(res);
}

maxabc <- function(mu1,mu2,sigma1,sigma2){
       
       dff=0;
       for(a in seq(0,6,0.1)){
        for(b  in seq(-1,2,0.1)){
          for(c in seq(0,0.4,0.05)){
              newdiff=pmean_diff(a,b,c)
              print(c(a,b,c,newdiff));
              if(newdiff > dff){dff=newdiff;aa=a;bb=b;cc=c;maxdiff=newdiff}
            }
        }
      }
      #maxLik(pmean_diff,start=c(0.2,1.2,0.1),constraints=list(param[3]<0.4,param[2]<2))
      
}
